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1.    Introduction 

In  these  lectures  we  will  survey  adaptive  methods  for  solving  partial  differential  equa- 
tions related  to  fluid  mechanics.  The  reader  will  find  however  that  we  spend  most  of  the 
time  discussing  issues  related  to  our  own  approach  -  adaptive  mesh  refinement.  Adaptive 
methods,  particularly  our  own,  are  rather  more  complicated  than  non-adaptive  methods.  Our 
approach  requires  internal  boundary  conditions  that  must  be  conservative,  data  structures  for 
keeping  track  of  several  layers  of  fine  grid  patches,  error  estimation,  and  heuristics  for 
automatic  grid  generation.  However,  all  this  machinery  fits  into  about  3,000  lines  of  FOR- 
TRAN code.  The  inner  loops  still  vectorize  and  dominate  the  total  execution  time  (our  own 
code  was  clocked  at  90  Megaflops  on  the  Cray  XMP).  In  practical  calculations  we  achieve 
gains  in  efficiency  by  factors  up  to  10  over  non-adaptive  methods.  By  citing  these  numbers, 
we  hope  to  reassure  the  intimidated  reader  and  disarm  the  skeptic.  Sophisticated  adaptive 
techniques  can  be  very  practical  tools  in  solving  the  most  difficult  flow  problems. 

The  aim  of  adaptive  calculations  is  to  put  greater  resolution  (i.e.  more  grid  points)  into 
regions  that  need  them.  By  now  it  is  accepted  that  this  is  a  useful,  if  not  essential,  thing  to 
do.  Often,  a  flow  will  be  generally  smooth  with  large  variations  (shocks,  boundary  layers, 
etc.)  occurring  in  small  subsets  of  the  domain.  It  is  usually  out  of  the  question  to  discretize 
an  entire  domain  using  the  finest  mesh  spacing  that  will  be  needed  throughout  an  entire  cal- 
culation. Indeed,  much  of  the  work  in  grid  generation  has  been  motivated  by  the  desire  to 
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cluster  grid  points  in  specific  places,  where  greater  accuracy  is  thought  to  be  needed. 

The  next  logical  step  is  to  automate  the  generation  of  an  adapted  grid.  The  finest  mesh 
spacing  needed  may  not  be  easy  to  guess  a  priori.  For  nonlinear  problems,  even  the  regions 
needing  extra  resolution  (e.g.  shocks)  may  not  be  known  in  advance.  An  automatic  adaptive 
method  is  needed  to  decide  where  to  put  the  points  during  the  computation  itself,  by  looking 
at  the  solution. 

We  can  distinguish  two  essentially  distinct  kinds  of  adaptive  methods.  The  first  method 
maintains  a  logically  rectangular  mesh,  but  deforms  the  mesh  so  that  grid  points  are  concen- 
trated in  regions  of  interest.  We  call  these  methods  Moving  Grid  Point  methods.  A  second 
approach  is  to  refine  computational  cells  by  adding  extra  grid  points  to  cells  that  lack  suffi- 
cient resolution.  Within  this  category  our  own  approach  is  to  organize  the  refined  cells  into 
logically  rectangular  fine  grid  patches  that  cover  regions  of  high  error.  We  call  this  the 
Embedded  Grid  approach.  We  now  discuss  these  methods  in  more  detail.  Most  of  our  dis- 
cussion of  these  approaches  assumes  they  are  used  in  conjunction  with  a  finite  difference 
scheme  for  solving  partial  differential  equations  (pdes). 

Moving  grid  point  methods  try  to  get  the  most  accurate  solution  for  a  fixed  cost.  They 
do  this  by  trying  to  optimally  distribute  a  fixed  number  of  grid  points  in  a  domain.  We  can 
broadly  recognize  two  methods  for  distributing  the  grid  points,  although  the  methods  can  cer- 
tainly be  combined.  In  fact,  some  of  the  most  effective  computations  have  used  combinations 
of  these  techniques.  In  one  approach,  grid  points  are  moved  at  each  time  step  by  integrating 
an  extra  system  of  equations  describing  the  grid  point  motion.  This  can  either  be  coupled 
with  the  original  system,  as  in  the  Moving  Finite  Element  method,  or  it  can  be  separate,  as  in 
a  Lagrangian  method.  An  advantage  of  having  grid  points  move  with  time  is  that  in  the  mov- 
ing coordinate  system,  the  solution  may  not  be  as  rapidly  changing  as  it  is  in  the  stationary 
coordinate  system.  Thus,  larger  time  steps  are  possible  without  sacrificing  accuracy.  How- 
ever, [Coyle,  Flaherty  and  Ludwig]  have  shown  that  some  reasonable  looking  equations  for 


mesh  movement  are  in  fact  unstable,  and  lead  to  mesh  trajectories  that  oscillation  badly. 

In  another  approach,  grid  points  are  redistributed  "between"  time  steps  according  to  a 
given  weight  function.  A  theoretical  motivation  for  this  approach  is  in  [Pereyra  and  Sewell], 
which  proves  the  optimality  in  one  dimension  of  equidistributing  a  mesh  based  on  the  trunca- 
tion error  of  the  solution.  In  practise,  however,  the  weight  functions  most  commonly  used 
are  spatial  derivatives  of  the  solution,  not  necessarily  related  to  the  truncation  error.  In  addi- 
tion, it  is  not  clear  how  to  generalize  this  to  many  dimensions. 

Moving  grid  point  methods  in  both  categories  often  have  trouble  maintaining  a  smooth 
grid.  Regularity  terms  and  penalty  functions  used  to  regularize  the  grid  add  overhead  and 
reduce  the  simplicity  and  robustness  of  these  methods.  This  problem  has  impeded  the  exten- 
sion of  moving  grid  point  methods  to  three  dimensions. 

We  will  describe  some  specific  methods  in  more  detail.  This  is  not  an  exhaustive  list  by 
any  means.  The  examples  are  chosen  to  illustrate  the  different  approaches. 

The  Moving  Finite  Element  (MFE)  method,  a  generalization  of  the  classical  finite  ele- 
ment method,  is  developed  in  [Miller  and  Miller;  Gelinas,  Doss,  Miller],  Theoretical  justifi- 
cation for  the  method  was  provided  by  [Dupont].  The  MFE  allows  the  positions  of  the  nodes 
defining  the  basis  elements,  s,,  as  well  as  the  element  amplitudes  a  ,  to  vary  in  time.  The 
idea  is  to  choose  the  s,  and  a,  to  minimize  (an  estimate  of)  the  error.  An  approximate  solu- 
tion v  to  u,  =  L(u)  depends  on  In  time  dependent  parameters 

"  =  v(°i(') *,('). *,('),  •  •  •  s„(t)), 

so  that 

v  -   - — a,  +    •  •  ■  - — a„  +   - — sx  +     ■  ■  ■  -2— , 
da  i  dan  dSl  ds„ 


=  £a>,  +  ifa. 


where 

3v  dv 

a .      3,   -    . 

da;  ds. 

The  norm  squared  of  the  residual  R   =   v-L(v)  is 

(R,R)  =   (v.w)  -  2(v,L(v))  -   (t(v),L(»)). 
This  expression  is  minimized  over  a  and  i  when 

—  (R,R)  =    —(R.R)  =  0,     i  =   1 n. 

This  gives  2n  ODE's  for  a,,  and  s  ,  of  the  form 

2 («..«;)«;  -  2(«,  P,)j.  -  (o„l(v))  =  0 

2(S„a,)fl;  +  2(3, ,(3,)  /,.   -   (p„L(v))  =  0, 

;  j 

and  so  the  solution  to  the  ODE's  tells  how  to  move  the  nodes  and  advance  the  solution. 

The  idea  behind  the  MFE  is  very  clever,  and  several  impressive  computations  have  been 
done  with  it.  But  as  with  other  moving  grid  point  methods,  to  keep  nodes  from  tangling, 
crossing  one  another,  or  bunching  in  one  region  leaving  holes  in  another  region,  various 
penalty  functions  or  ad  hoc  regularizing  terms  have  to  be  added.  These  terms  so  far  are  not 
very  satisfying.  In  addition,  the  system  of  ODE's  is  a  non-linear  implicit  system  which  is 
often  very  stiff.  Much  of  this  stiffness  is  due  to  the  s  terms:  the  adaptivity  has  made  the 
equations  harder  to  solve. 

A  one  dimensional  mesh  is  said  to  be  equidistributed  with  respect  to  a  weight  function  w 
if  w,  ~  constant  at  every  grid  point.  Typical  weight  functions  used  involve  the  gradient  of 
the  solution,  a  second  derivative,  or  combinations  of  these.  Examples  of  this  approach  are 
found  in  [Dwyer,  Kee  and  Sanders;  Dwyer].  In  their  work,  the  grid  points  are  distributed 
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along  one  coordinate  direction  so  that 

/(I   -0,1^1   -b2\^\)ds 
J0  ds  ds- 

7,  ds  ax- 

is constant.  Here,  7"  is  a  dependent  variable  that  reflects  the  variation  in  the  solution 
(several  can  be  used),  s  is  arclength,  and  the  constants  b\,bi  are  picked  by  the  user.  For  two 
dimensional  problems,  there  is  hopefully  a  clear  coordinate  direction  along  which  to  adapt 
the  grid.  For  many  problems,  this  technique  works  well,  although  the  authors  report  some 
problems  with  grid  oscillations  and  grid  skewness. 

A  different  mesh  equidistribution  criterion  has  been  used  in  one  dimension  by  [Winkler] 
to  calculate  spherically  symmetric  star  formation.  Since  the  dependent  variables  vary  up  to 
20  orders  of  magnitude,  Winkler  uses  logarithms  in  his  weight  function, 

2  f»  (ln-(^-)-ln3(^-))  =  0. 

h  in  \  anabies  .'  J 

to  bt  ttjuijisfiDured 

Equidistribution  of  the  dependent  variables  implicitly  defines  the  mesh  motion.  These  equa- 
tions are  then  coupled  with  the  pdes  for  the  original  time  dependent  variables.  Everything  is 
advanced  in  time  using  an  implicit  method  and  Newton  iteration. 

However,  in  more  recent  work  [Winkler,  Mihalas,  Norman]  the  grid  motion  is  more 
explicitly  accomplished  by  solving  an  elliptic  equation  for  the  grid  motion 

-^  =  (A/r),  -  (A/5).. 
of 

Here,  there  are  two  competing  forcing  functions:  fr  controls  grid  smoothness,  and  fs  controls 
the  resolution  in  the  solution.  The  actual  definition  of  these  functions  is  quite  complex  and 
rather  ad  hoc.  However,  with  proper  choice  of  numerical  parameters,  the  method  can  per- 
form very  well. 


A  different  approach  to  equidistribution  is  developed  by  [Rai  and  Anderson;  Anderson 
and  Rai],  They  try  to  directly  compute  grid  point  speeds,  and  then  integrate  to  determine  the 
coordinates  at  the  next  time  step.  The  goal  is  to  equidistribute  a  positive  quantity  e  over  the 
mesh,  so  that  e,  ~  constant  ~  ea,  at  each  grid  point.  This  quantity  e  is  hopefully  correlated 
with  the  error.  Regions  where  et>ea  need  more  points,  and  vice  versa.  This  is  done  by  hav- 
ing a  grid  point  x,  attract  or  repel  a  grid  point  x,  by  inducing  a  velocity  proportional  to 
e,  —  eav,  and  inversely  proportional  to  the  distance  |x, -x .  |,  on  the  grid  point  x,.  To  get  the 
velocity  at  one  point,  the  contributions  of  the  other  n  -  1  grid  points  are  summed.  The  idea  is 
extended  to  rwo  dimensions  by  determining  a  grid  point  speed  in  each  coordinate  direction 
separately.  Boundary  points  are  constrained  to  stay  on  the  boundary  by  only  allowing  motion 
in  the  direction  tangential  to  the  boundary.  Grid  velocities  are  kept  below  a  user  specified 
maximum  by  scaling  the  constant  of  proportionality  at  every  step.  This  maximum  values  must 
be  chosen  carefully  -  large  values  lead  to  grid  oscillations  and  small  values  inhibit  grid  move- 
ment. Additional  limits  are  imposed  by  various  damping  techniques  to  prevent  excessive 
mesh  stretching  and  grid  point  clustering.  On  test  problems  with  a  single  discontinuity  this 
method  works  well.  Without  more  grid  control,  it  is  not  clear  how  well  this  method  would 
perform  on  complicated  shock  interaction  problems,  where  highly  skewed  meshes  could 
result. 

A  nicely  formulated  and  more  systematic  method  for  adapting  a  grid  is  by  [Brackbill 
and  Saltzman;  Brackbill;  Saltzman  and  Brackbill].  They  define  three  measurable  properties  of 
a  computational  mesh:  its  smoothness,  orthogonality,  and  cell  volume  variation.  This  latter  is 
used  to  measure  the  accuracy  of  the  grid;  cell  volumes  should  be  small  when  a  user  specified 
weight  function  such  as  the  solution  gradient  is  large.  An  optimal  mesh  optimizes  these  three 
properties.  For  example,  if  £,-n  represent  the  coordinates  of  the  computational  mesh  in  two 
dimensions,  smoothness  can  be  measured  by 

ls=   $  (?£):   -    (Vt,)2  dV, 
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orthogonality  by 

/o  =  J  (Vt-Vil)2 /3  dV, 

D 

and  volume  variation  by 

/    -  /  w(x,y)  J  dV, 

D 

where  w  is  user  specified  ,  and  J  is  the  Jacobian  of  the  mapping  J  =    — —      -   - —  ,,       The 

adaptive  mesh  is  generated  by  optimizing  a  linear  combination 

/  =  Is  +  M,  *  XB/0, 

where  \,  ,\,,  are  positive  constants.  Finding  this  mesh  is  a  variational  problem  for  £  and  i\. 
For  time  dependent  pdes,  a  few  Gauss  Seidel  iterations  between  time  steps  suffices  to  define 
the  new  mesh  at  the  next  time  step.  The  solution  can  then  be  interpolated  (conservatively  if 
necessary)  onto  the  new  mesh. 

As  with  all  methods,  some  knobs  must  be  tweaked  for  proper  performance.  Here,  for 
example,  the  authors  report  some  smoothing  must  be  performed  on  the  weight  function  w 
itself,  the  maximum  and  minimum  cell  sizes  must  be  constrained,  and  the  /3  is  needed  in  the 
orthogonality  integral  above  to  reduce  rounding  error  problems. 

Several  people  are  applying  this  technique  to  fluid  flow  problems.  [Kreis,  Thames  and 
Hassan]  use  this  variational  approach  to  compute  steady  two  dimensional  flow  around  an  air- 
foil. [Nakahashi  and  Deiwert]  apply  a  simplified  but  somewhat  more  ad  hoc  version  of  this 
approach  in  solving  complicated  viscous  and  inviscid  flow  problems  in  two  and  three  dimen- 
sions. [Yanenko,  et  alj  have  also  used  variational  methods  for  moving  meshes  to  solve  prob- 
lems in  gas  dynamics. 

One  drawback  of  any  moving  grid  point  method  is  the  time  step  restriction.  When  using 
an  explicit  finite  difference  scheme  to  integrate  the  solution,  the  smallest  mesh  spacing  any- 
where in  the  grid  determines  the  allowable  time  step  for  the  entire  grid.    This  global  penalty 


causes  some  users  of  moving  grid  point  methods  to  set  a  minimum  cell  size,  although  this 
reduces  the  adaptivity.  In  general,  however,  moving  grid  point  methods  have  the  potential  to 
achieve  a  much  finer  grid  spacing  with  fewer  points  than  the  more  structured  embedded  grid 
methods  can.  Finally,  the  effect  of  very  rapid  grid  variations  and  grid  skewness  on  solution 
accuracy  is  not  well  understood. 

In  the  embedded  grid  approach  to  adapting  a  grid,  greater  resolution  is  achieved  by 
adding  new  grid  points,  not  by  clustering  existing  grid  points.  In  moving  grid  point  methods, 
if  there  are  several  "hot  spots"  in  the  solution,  the  grid  points  tend  to  cluster  in  the  worst 
region,  leaving  the  other  spots  unresolved.  With  local  grid  refinement,  if  there  are  several 
areas  that  need  more  resolution,  points  are  added  to  all  of  them.  Whereas  moving  grid  point 
methods  try  for  the  best  solution  at  a  fixed  cost,  local  grid  refinement  methods  try  to  attain  a 
fixed  accuracy  for  a  minimum  cost. 

We  can  distinguish  two  approaches  to  local  grid  refinement.  In  one  approach,  grid 
points  are  added  by  refining  existing,  coarser  grid  cells  on  a  cell  by  cell  basis.  This  approach 
is  commonly  found  in  adaptive  finite  element  calculations.  In  the  other  approach,  the  addi- 
tional grid  points  are  organized  into  fine  grid  patches  that  are  overlaid  on  the  coarse  grid. 
This  allows  the  possibility  of  rotating  the  fine  grids  with  respect  to  the  underlying  coarse 
grid,  so  that  the  grid  is  aligned  with  some  feature  in  the  solution.  Figure  1.1  illustrates  a  typi- 
cal refinement  pattern  of  both  approaches. 

In  the  first  approach,  it  is  possible  to  refine  only  those  coarse  cells  needing  refinement. 
However,  since  this  usually  forms  an  irregular  pattern,  linked  lists  and  other  data  structures 
are  used  so  that  each  cell  knows  if  the  neighboring  cell  is  refined  or  not.  This  slows  down 
the  integrator  and  makes  vectorization  more  difficult.  In  the  second  approach,  the  cells  need- 
ing refinement  are  organized  into  rectangular  grid  patches.  This  means  that  cells  that  don't 
need  more  resolution  are  refined  along  with  the  cells  that  do.  Heuristic  grid  generation  algo- 
rithms  are   needed    to   create    the    fine    grid    patches.    However,    since    the    fine    grids    are 
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Fig.  1.1:  in  Fig.  1.1a  individual  coarse  grid  cells  are  refined;  in  Fig.  1.1b  entire  fine  grid 
patches  are  superimposed  on  the  coarse  grid. 


themselves  logically  rectangular,  an  off-the-shelf,  vectorizing  integrator  can  be  used  for  both 
the  coarse  and  fine  grids. 

One  benefit  of  either  approach  to  local  grid  refinement  is  that  only  the  finest  grid  cells 
need  to  have  the  smallest  time  step.  Unrefined  grid  cells  can  use  the  original  time  step. 
Local  grid  refinement  methods  thus  refine  the  grid  in  both  space  and  time.  For  hyperbolic 
pdes  it  is  common  to  refine  the  space  step  and  time  step  by  the  same  factor.  The  advantage 
of  this  is  that  an  explicit  finite  difference  scheme  is  stable  regardless  of  whether  a  cell  has 
been  refined,  since  the  mesh  ratio  of  time  step  to  space  step  is  constant  for  all  cells.  A  disad- 
vantage of  local  grid  refinement  methods  is  that  special  difference  equations  must  be  used  at 
the  interfaces  between  refined  and  unrefined  cells. 

Embedded  grid  methods  are  usually  "recursive"  -  if  a  refined  region  still  has  a  high 
error,   it  is  refined  again.    This  leads  to  a  natural  grid  hierarchy,  of  coarse  to  fine  grids. 
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Where  embedded  grid  methods  in  both  categories  differ  is  in  how  they  structurally  organize 
the  refined  regions.  In  addition,  the  criteria  used  by  embedded  grid  practitioners  to  decide 
where  to  refine  the  grid  is  as  varied  as  it  is  in  moving  grid  points  methods.  To  a  large  extent 
however,  the  differences  within  each  category  are  not  fundamental,  but  are  due  to  the 
amount  of  effort  spent  implementing  an  algorithm.  Embedded  grid  refinement  is  also 
becoming  popular  in  multigrid  calculations  [Bai  and  Brandt;  McCormick  and  Thomas], 
although  these  are  not  yet  adaptive. 

We  will  briefly  describe  several  approaches  to  local  grid  refinement  in  both  categories. 
The  rest  of  the  notes  are  spent  describing  our  own  approach  to  local  grid  refinement  using 
fine  grid  patches.  We  note  that  many  of  the  original  ideas  for  local  grid  refinement  for 
hyperbolic  equations  were  first  discussed  in  [Oliger].  A  stability  theory  for  local  grid  refine- 
ment for  hyperbolic  equations  is  developed  in  [Ciment(1972);  Berger(1985)]. 

Preliminary  work  demonstrating  the  feasibility  of  local  mesh  refinement  was  done  by 
[Bolstad]  in  one  dimension,  and  [Gropp(1981)]  in  two  dimensions.  In  Gropp's  implementa- 
tion, the  computational  cells  of  a  coarse  grid  needing  refinement  were  organized  into 
columns,  from  a  starting  row  to  an  ending  row.  Two  disjoint  segments  of  the  same  column 
could  not  be  refined  without  refining  cells  in  between.  Also,  he  did  not  allow  recursive 
refinement  of  cells.  More  recent  work  is  recursive.  [Gropp(1986)|  uses  cell  by  cell  refine- 
ment while  maintaining  a  "1-regular"  grid,  where  no  more  than  one  irregular  grid  point  is 
allowed  per  cell  side.  This  means  that  if  a  coarse  cell  is  refined  twice,  the  neighboring  coarse 
cell  must  also  be  refined  at  least  once.  This  is  done  to  insure  that  the  work  at  the  interface 
between  refined  and  unrefined  cells  doesn't  get  too  complicated,  at  the  price  of  refining  more 
cells  than  absolutely  necessary.  Also,  the  mesh  refinement  factor  here  is  limited  to  two. 

Cell  by  cell  refinement  has  been  used  by  [Usab  and  Murman]  in  their  simulations  of 
flow  past  an  airfoil.  Here  too  the  mesh  refinement  factor  is  two,  but  since  the  integrator  uses 
a  multigrid  strategy  to  accelerate  convergence  to  a  steady  state  solution,  this  is    a  natural 
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refinement  ratio.  The  data  structure  for  this  type  of  refinement  uses  V  pointers  per  cell,  but 
has  no  topological  restriction  on  the  shape  of  the  domain  or  the  refined  cells.  Also,  this 
refinement  is  not  automatic;  the  refined  cells  must  be  established  at  the  start  of  a  computa- 
tion. 

A  similar  adaptive  strategy  was  developed  by  [Dannenhotfer  and  Baron).  They  experi- 
ment with  using  first  and  second  differences  of  density,  pressure,  velocity  and  entropy  to 
select  the  cells  needing  refinement.  The  first  difference  of  density  works  best  here.  An 
interesting  development  in  this  work  is  a  method  to  automatically  choose  the  threshold  value 
above  which  cells  are  refined.  By  considering  the  percentage  of  refined  cells  as  a  function  of 
refinement  tolerance,  a  threshold  value  is  chosen  above  the  value  at  which  half  the  points  are 
refined,  but  hopefully  below  the  value  needed  to  distinguish  special  features  in  the  solution. 

Our  own  work  on  adaptive  grid  refinement  uses  fine  grid  patches  in  regions  requiring 
greater  resolution.  These  patches  can  be  rotated  [Berger  and  Oliger],  or  aligned  with  and 
anchored  at  coarse  grid  grid  points  [Berger  and  Jameson;  Berger  and  Colella].  Aligned  grids 
are  used  for  problems  with  shocks,  since  for  these  problems,  the  difference  equations  at  the 
interface  between  fine  and  coarse  grids  should  be  conservative.  This  is  easier  to  do  if  the  fine 
and  coarse  grids  are  in  the  same  coordinate  system.  We  use  Richardson-type  estimates  of  the 
local  truncation  error  to  find  the  regions  which  need  to  be  refined.  A  tree  data  structure  is 
used  to  keep  track  of  the  grids.    A  complete  presentation  of  this  approach  starts  in  section  2. 

For  the  hardest  computational  problems,  the  best  approaches  will  undoubtedly  use  com- 
binations of  these  techniques.  By  allowing  grid  points  to  move,  a  mesh  can  be  aligned  with 
some  feature  in  the  solution  without  the  extra  overhead  of  using  rotated  fine  grid  patches  to 
do  the  alignment.  The  use  of  local  grid  refinement  to  meet  an  accuracy  criterion  can  prevent 
the  global  time  step  restrictions  of  moving  grid  point  methods. 
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Some  work  in  this  direction  has  already  begun.  In  [Gropp  (1984)|,  a  method  of  local 
grid  refinement  using  fine  grid  patches  was  developed  for  solving  hyperbolic  pdes.  Each 
patch  is  allowed  to  move  at  a  uniform  speed  that  changes  at  each  time  step.  In  test  problems 
to  date,  this  speed  must  be  set  by  the  user,  or  calculated  as  part  of  the  solution  by  averaging 
the  velocity  over  all  grid  points  in  a  fine  grid. 

In  [Adjerid  and  Flaherty),  local  refinement  was  combined  with  a  moving  finite  element 
method  for  parabolic  pdes  in  one  space  dimension.  The  refinement  algorithm  adds  elements 
to  subintervals  based  on  local  error  estimates.  In  addition,  all  mesh  lines  move  according  to  a 
mesh  equidistribution  strategy.  A  variety  of  techniques  are  combined  in  [Arney  and  Flaherty) 
in  extending  this  work  to  solve  hyperbolic  problems  in  two  space  dimensions.  Here,  cell  by 
cell  refinement  is  coupled  with  a  moving  base  mesh.  They  clearly  show  that  moving  meshes 
can  improve  the  efficiency  of  local  grid  refinement. 

2.    The  Adaptive  Mesh  Refinement  Algorithm 

We  have  two  reasons  for  preferring  to  use  locally  uniform  fine  grid  patches.  The  first 
is  simplicity.  Rectangular  grids  have  the  simplest  user  interface.  We  use  the  same  integrator 
on  the  fine  grid  as  on  the  coarse  grid.  By  separating  the  integrator  from  the  adaptive  stra- 
tegy, an  off-the-shelf  integrator  can  be  used  without  modification.  This  eliminates  much  of 
the  problem  specific  work  in  doing  adaptive  calculations  (see,  however,  the  section  on  inter- 
face conditions).  The  uniform  structure  also  leads  to  simple  vectorization  and  parallelization 
of  adaptive  codes.  The  grid  management  (see  the  section  on  data  structures  below)  overhead 
in  our  approach  is  proportional  to  the  number  of  grid  patches,  (on  the  order  of  10,  say), 
rather  than  the  number  of  refined  cells. 

The  second  reason  for  favoring  locally  uniform  patches  is  that  wave  propagation  proper- 
ties are  much  different  on  uniform  than  on  non-uniform  meshes.  For  example,  [Giles  and 
Thompkins)  have  shown  that  on  highly  stretched  grids  around  airfoils,  waves  can  be  trapped 


•  13  - 

close  to  the  body,  rather  than  propagating  towards  the  far  field  as  they  should.  Stretched 
grids  can  have  a  particularly  adverse  effect  on  shock  wave  propagation  [P.  Colella,  personal 
communication). 

Let  us  now  discuss  the  specific  pieces  of  our  approach  as  they  stand  at  present.  These 
will  be  examined  in  detail  below.  We  use  Richardson  extrapolation,  a  procedure  for  estimat- 
ing truncation  error,  to  decide  where  refinement  is  needed.  Then  we  apply  grid  generation 
heuristics  to  create  fine  grid  patches  covering  the  regions  that  need  refinement.  Some  of 
these  heuristics  are  related  to  clustering  strategies  used  in  computer  vision  and  artificial  intel- 
ligence. A  tree  data  structure  stores  the  information  describing  the  nested  grids,  allowing 
access  to  the  individual  grid  patches  as  needed.  Finally,  we  must  develop  numerical  boun- 
dary conditions  for  the  interfaces  between  fine  and  coarse  grid  patches.  Unfortunately,  it  is 
not  always  easy  to  find  boundary  conditions  that  are  stable,  conservative,  and  minimize 
unwanted  wave  reflection. 

We  have  been  careful  to  make  our  adaptive  mesh  refinement  algorithm  as  automatic  as 
possible.  The  user  must  supply  the  integration  subroutine,  the  coarse  grid  boundary  condi- 
tions, and  a  subroutine  describing  the  domain.  The  only  essential  input  parameter  the  user 
must  specify  is  the  error  tolerance  that  triggers  the  grid  refinement. 

When  all  these  components  are  assembled,  a  typical  integration  step  proceeds  as  fol- 
lows. The  integration  steps  on  different  grids  are  interleaved  so  that  before  advancing  a  grid, 
all  the  finer  level  grids  have  been  integrated  to  the  same  point.  Before  every  coarse  grid 
step,  all  grids  have  been  integrated  to  the  same  time.  One  coarse  grid  cycle  is  then  the  basic 
unit  of  the  algorithm.  Below,  we  outline  a  coarse  grid  integration  cycle  as  a  recursive  pro- 
cedure. Of  course,  in  writing  the  mesh  refinement  program  in  FORTRAN,  we  have  con- 
verted it  to  an  iterative  procedure.  The  variable  r  denotes  the  mesh  refinement  factor  in  time. 
The  regridding  procedure  is  done  every  few  steps,  so  any  particular  step  may  or  may  not 
involve  regridding. 
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Recursive  Procedure  Integrate  (level) 
Repeat  r'"'1  times 

Regridding  time  ?  -     error  estimation  and  grid  generation  for  grids  at  level  level   and  finer 
step  Af/Plf/  on  all  grids  at  level  level 
if  (level+  1)  grids  exist 
then  begin 

integrate  (level-  1) 
updatc(tevel,  level— \) 
end 
end. 

level  =  0    (*  coarsest  grid  level  *) 
Integrate  (level) 

Coarse  Grid  Integration  Cycle 

The  update  procedure  above  needs  further  explanation.  If  a  fine  grid  is  nested  in  a 
coarser  grid,  then  when  they  are  integrated  to  the  same  time,  the  coarse  grid  values  are 
updated  by  injecting  the  fine  grid  solution  values  onto  the  coarse  grid  points.  If  grid  points 
do  not  match  up  at  regular  intervals,  for  example,  if  the  grids  are  rotated  with  respect  to  each 
other,  interpolation  is  used  to  find  the  value  from  the  fine  grid  that  replaces  the  coarse  grid 
point's  solution.  If  the  calculation  is  to  be  conservative,  then  we  do  not  (as  yet)  use  rotated 
rectangles,  and  the  updating  is  done  by  using  a  volume  weighted  average  of  all  fine  grid 
values  in  the  coarse  cell.  For  example,  if  the  indices  i,j  refer  to  the  coarse  cell  being 
updated,  and  the  indices  k,m  refer  to  the  lower  left  fine  grid  value  in  that  cell,  the  solution  is 
updated  by 


15  - 


r      p  «0  4=g 


The  updating  step  also  takes  care  of  conservatively  matching  the  difference  equations 
used  on  the  fine  and  coarse  grids  across  the  grid  interface.  This  will  be  described  in  the  sec- 
tion on  interface  conditions. 

3.    Error  Estimation 

The  decision  to  refine  the  grid  is  based  on  estimates  of  the  local  truncation  error,  as  was 
proposed  by  [Oliger].  This  is  motivated  by  convergence  results  of  [Gustafsson]  for  the  initial 
boundary  value  problem.  We  use  a  method  based  on  Richardson  extrapolation.  For  simpli- 
city, let  Q  be  a  two-level  explicit  finite  difference  operator,  h  the  mesh  width  in  space,  and  k 
the  time  step.  Suppose  that  we  know  that  the  method  we  are  using  has  accuracy  of  order  q  in 
both  time  and  space.  If  the  solution  is  smooth  enough,  the  local  truncation  error  is  defined 
as 

u{x,t-k)-Qu(x,t)  =  k  (  cx{x,t)kl  -  c2(x,t)hq)  -  *  0(kq]   -  h"']) 

=  t(x,0  +  k  0{kq^  -  h*-]), 

where  the  leading  term  is  denoted  by  t.  If  u  is  smooth  enough,  then  if  we  take  two  time  steps 
with  the  method  Q,  to  leading  order  the  error  is  2t  , 

u(x,t  +  2k)-Q2u(x,t)  =  2t  +  k  0(k"'1   *  hq~l) 

Let  Q2h  denote  the  same  difference  method  as  Q  but  based  on  a  mesh  widths  of  2/i  and  2k. 
Then 

u(x,r  +  2*)  -  Q2,u{x,t)  =  2"-'t  +  0(h«-z). 

By  taking  two  steps  with  the  regular  integration  scheme,  and  one  "giant"  step  using  every 
other  grid  point,  the  difference 

Q2u(x,t)  -  Q^,,u{x,t) 

, =     T    *    0(h°      -) 

2"  " '  -  2 
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can  be  formed,  giving  an  estimate  of  the  local  truncation  error  at  time  r  This  is  illustrated 
schematically  in  Figure  3.1  for  a  one  dimensional  difference  scheme  with  a  three  point  sten- 
cil. 


Fig.  3.1  Richardson  error  estimation  procedure  in  one  space  dimension  and  time. 


This  procedure  has  several  points  to  recommend  it.  First,  it  is  not  necessary  to  know  the 
exact  form  of  the  truncation  error.  The  c,  and  c2  above  involve  higher  derivatives  of  the 
solution,  but  need  not  be  known  explicitly.  Especially  for  higher  order  accurate  methods  and 
for  systems  of  equations  in  several  variables,  it  can  be  quite  difficult  to  compute  the  exact 
form  of  the  truncation  error.  Not  only  is  our  estimator  independent  of  the  pde,  the  error  esti- 
mation procedure  is  independent  of  the  difference  method,  as  long  as  the  accuracy  in  time 
and  space  are  the  same.  This  is  important  if  mesh  refinement  is  going  to  be  applicable  to  a 
wide  variety  of  problems  without  much  re-programming.  If  the  spatial  and  temporal  order  of 
accuracy  are  different,  then  a  slightly  more  complicated  variation  can  be  used. 
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For  problems  with  smooth  solutions,  the  Richardson  procedure  provides  a  robust  esti- 
mate of  the  local  truncation  error.  For  non-smooth  solutions,  we  no  longer  have  an  accurate 
error  estimate.  However,  the  Richardson  estimates  still  provide  a  good  criterion  for  refine- 
ment, since  near  a  singularity,  the  estimate  will  be  large.  For  piecewise  constant  initial  condi- 
tions, the  Richardson  algorithm  gives  an  estimate  proportional  to  the  jump  in  the  solution. 

4.    Grid  Generation 

If  there  are  several  nested  levels  of  refined  grids,  the  error  estimation  and  grid  genera- 
tion procedures  are  applied  on  each  level  to  (re-)create  the  fine  grids  at  the  next  level.  The 
error  estimation  procedure  produces  a  list  of  coarse  grid  points  with  high  error  estimates, 
indicating  that  a  fine  grid  patch  is  needed  in  that  region.  Our  grid  generation  algorithms  try 
to  produce  grids  that  have  as  little  overlap  as  possible,  so  that  the  area  that  is  unnecessarily 
refined  is  as  small  as  possible.  The  algorithm  also  strives  for  a  small  number  of  patches,  to 
reduce  overhead  and  overlap.  It  is  difficult  to  find  a  foolproof  algorithm  that  satisfies  these 
often  conflicting  goals.  Our  approach  is  to  start  with  a  simple  algorithm  that  works  in  most 
cases,  and  try  to  detect  when  it  does  a  bad  job  of  fitting  grids;  in  these  cases  we  use  a  more 
expensive  algorithm.    More  details  about  these  algorithms  can  be  found  in  [Berger  (1986)]. 

The  first  approach  we  use  clusters  the  grid  points  using  a  nearest  neighbor  algorithm. 
Each  cluster  will  later  become  a  fine  grid.  The  algorithm  starts  with  one  flagged  point  in  a 
new  cluster.  Successive  points  are  included  in  this  cluster  if  the  distance  from  the  point  to  the 
cluster  is  less  than  two  mesh  widths. 

The  nearest  neighbor  algorithm  is  very  good  at  separating  spatially  distinct  phenomena, 
so  that  different  features  of  the  solution  will  be  in  separate  grids.  For  example,  in  flow 
around  an  airfoil,  the  leading  and  trailing  edges  of  the  airfoil  will  typically  need  refinement. 
The  flagged  coarse  grid  cells  with  high  error  estimates  will  naturally  separate  into  two  clus- 
ters. However,  the  algorithm  fails  when  there  is  one  long  feature  needing  refinement,  such  as 
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a  bow  shock,  that  is  not  aligned  with  a  grid  line  of  the  coarse  grid.  In  these  cases,  unless  a 
rotated  rectangle  is  used,  one  large  refined  unrotated  grid  patch  containing  all  the  flagged 
points  would  include  too  many  unnecessarily  refined  cells  as  well,  and  would  be  inefficient. 
For  more  complicated  shock  interaction  patterns,  found  for  example  in  mach  stem  computa- 
tions, one  rotated  rectangle  would  still  be  inefficient.  Several  contiguous  fine  grid  patches 
are  needed,  as  illustrated  in  Figure  4.1.  This  inefficiency  is  detected  by  counting  the  number 
of  flagged  points  versus  the  total  number  of  coarse  cells  in  the  new  fine  grid.  If  this  ratio 
falls  below  50%,  a  more  expensive  clustering  algorithm  must  be  used. 


Fig.  4.1  One  large  grid  patch  would  refine  too  much  of  the  grid.  Several  grid  patches 
are  more  efficient. 


We  have  developed  two  alternative  approaches  for  the  more  complicated  clustering 
case.  The  starting  point  for  both  algorithms  is  the  inefficient  nearest  neighbor  cluster.  In  the 
first  algorithm,  we  use  data  structures  to  organize  the  flagged  points  within  the  cluster,  so 
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that  a  smart  subdivision  of  the  points  can  be  made.  The  flagged  points  are  connected  into  a 
minimal  spanning  tree  (MST).  An  MST  is  the  connected,  acyclic  graph  connecting  all  points 
so  that  the  sum  of  the  lengths  of  the  edges  is  a  minimum.  In  this  graph,  each  flagged  point  is 
connected  to  is  nearest  neighbor.  The  tree  must  now  be  disconnected  by  breaking  edges  so 
that  several  subclusters  are  formed.  This  can  be  done  iteratively.  Start  with  one  point  in  a 
new  subcluster.  Add  neighboring  points  in  the  MST  as  long  as  the  resulting  grid  is  still  effi- 
cient. This  is  the  bottom-up  approach  to  grid  generation. 

The  other  approach  is  top-down:  use  a  bisection/merging  algorithm  on  grid  patches  that 
are   inefficient.   The  long  direction   of  the   inefficient  rectangular  grid  is  bisected,   and  the 
flagged  points  are  sorted  into  two  clusters  depending  on  which  half  they  are  in.  The  process 
is  repeated  on  the  two  clusters.  The  bisection  step  ends  when  each  cluster  has  an  acceptable 
efficiency   raring.   However,  since  no  geometric   information   is  used   in  this  approach,   the 
bisection  step  can  accidentally  create  several  small  clusters  that  would  be  more  efficient  if 
joined  together.  For  this  reason,  the  bisection  step  is  followed  by  a  merge  step.  Two  potential 
fine  grids  are  merged  if  the  new  grid    would  be  more  efficient  or  cost  less  than  the  two 
smaller  grids.    The  cost  function  we  use  to  measure  this  is  proportional  to  the  cost  of  an 
integration  step  on  each  grid,  which  fully  vectorizes,  plus  the  cost  associated  with  the  perime- 
ter of  each  grid,  which  often  uses  scalar  arithmetic.    The  merging  step  ends  when  no  pair  of 
grids  can  be  successfully  merged.    Here,  too,  more  sophisticated  data  structures  could  play  a 
role,  by  only  attempting  to  merge  clusters  that  are  connected  in  an  MST  for  example.    This 
algorithm  was  used  to  create  the  grids  in  Figure  4.2. 

Given  a  cluster  of  flagged  points,  it  is  easy  to  determine  the  dimensions  of  the  fine  grid 
patch  so  that  it  contains  all  the  points.  If  a  rotated  rectangle  is  desired,  it  is  also  easy  to  deter- 
mine the  angle  of  rotation  of  the  grid  patch.  The  subgrid  should  have  the  same  orientation  as 
the  flagged  points.  We  approximate  this  using  the  first  and  second  moments  of  the  flagged 
points.  Specifically,  in  two  dimensions,  let  x,   y,    1  <  i  <  n  be  the  coordinates  of  the  flagged 
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Fig.   4.2  The  bisection/merge  algorithm   produced   the   fine  grid   patches  around  the 
discontinuity  shown  in  the  contour  plot. 


points,  and  x,  y  their  mean.  Consider  the  symmetric  matrix 


M'M  = 


£  *,y<  -  xy  S  y:  -  y" 


where 


M  = 


I  I 


x,  -  x    y,  -  y 


This  matrix  determines  an  ellipse  with  the  same  first  and  second  moments  as  the  flagged 
points  [Cramer].  The  axes  of  the  ellipse  are  the  eigenvectors  of  the  2  by  2  matrix  M'M, 
which  we  use  to  determine  the  orientation  of  the  rectangle.  Since  the  matrix  is  only  2  by  2, 
regardless  of  the  number  of  flagged  points,  it  is  easy  to  find  the  eigenvectors.    The  most 
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expensive  part  of  this  algorithm  is  the  determination  of  the  dimensions  of  the  rectangle,  given 
its  orientation,  so  that  all  flagged  points  are  included.  It  turns  out  this  algorithm  is 
equivalent  to  doing  a  total  least  squares  fit  to  the  flagged  points,  finding  the  slope  of  the  line 
through  the  mean  of  the  points  that  minimizes  the  sum  of  the  squares  of  the  distances  from 
each  point  to  the  line. 

To  finish  the  discussion  of  the  complete  grid  generation  algorithm,  several  annoying 
details  must  be  mentioned  here.  Care  must  be  taken  to  insure  that  fine  grids  are  properly 
nested  in  the  next  coarser  grids.  In  effect,  this  means  the  fine  grid  error  estimates  are  used 
instead  of  the  coarse  grid  estimates  at  the  same  point.  This  is  done  by  flagging  coarse  grid 
points  corresponding  to  already  created  two-level-finer  grid  patches,  thus  insuring  the 
existence  of  an  intermediate  level  fine  grid.  Secondly,  a  buffer  zone  of  unflagged  points  is 
added  around  every  grid.  This  insures  that  discontinuities  or  other  regions  of  high  error  do 
not  propagate  out  from  a  fine  grid  into  coarser  regions  before  the  next  regridding  time.  The 
larger  the  buffer  zone,  the  more  time  steps  can  be  taken  to  lengthen  between  regnddings. 
This  buffer  zone  is  added  by  flagging  all  coarse  grid  points  that  are  sufficiently  close  to 
flagged  points  with  high  error  estimates.  A  buffer  zone  of  2  cells  in  each  direction  is  typical. 
By  flagging  neighboring  points,  instead  of  enlarging  the  grids  at  a  later  step,  the  area  of  over- 
lap between  grids  is  reduced.  Third,  if  the  grids  at  a  given  level  do  not  form  a  convex 
region,  it  may  happen  that  finer  level  grids  stick  out.  This  is  an  accident  of  the  grid  genera- 
tion only.  Proper  nesting  of  grids  is  enforced  by  repeatedly  subdividing  the  fine  grid  until 
each  piece  does  fit  in  the  coarser  grids.  Since  the  flagged  points  originally  were  inside  the 
non-convex  coarse  grids,  the  new  grids  containing  the  flagged  points  must  eventually  lie 
inside  as  well. 

Although  these  procedures  are  still  experimental  and  under  development,  they  have 
already  been  used  on  several  different  types  of  problems.  In  addition  to  solving  the  Euler 
equations,   the  grid  generation  algorithms  have   been   used   in  solving  steady   state   Navier 
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Stokes  problems  [Caruso,  Ferziger  and  Oliger),  as  well  as  pdes  in  meteorology  [Oliger  and 
Wijngaart,  to  appear).  Although  not  optimal,  they  have  been  proven  successful  and  efficient 
in  the  automatic  generation  of  adaptive  subgrids. 

5.    Data  Structures 

The  data  structures  in  our  adaptive  mesh  refinement  strategy  turn  out  to  be  surprisingly 
simple,  although  they  are  crucial  to  the  feasibility  of  the  entire  approach.  The  data  structure 
we  use  to  manage  the  grids  can  be  thought  of  as  a  tree,  where  each  node  in  the  tree 
represents  a  grid.  A  parent  node  in  the  tree  represents  a  coarser,  parent  grid.  Subgrids  are 
considered  the  offspring  of  their  parent  grid.  The  data  structure  is  not  exactly  a  tree,  how- 
ever, since  a  subgrid  can  have  more  than  one  parent  grid.  Figure  5.1  shows  a  typical  grid 
hierarchy,  and  the  linked  list  that  represents  it. 
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Fig.  5.1  shows  a  typical  2D  data  structure. 
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Because  this  tree  structure  can  grow  or  shrink  dynamically,  some  form  of  dynamic 
storage  allocation  is  needed,  both  for  the  grid  information  in  each  node  and  the  solution 
storage  for  each  grid.  Since  FORTRAN  does  not  provide  such  a  facility,  the  storage  manage- 
ment must  be  explicitly  provided  by  the  mesh  refinement  program.  We  keep  a  linked  list  of 
free  nodes  (up  to  a  maximum  of  100)  which  are  assigned  to  newly  created  grids,  and 
reclaimed  when  a  grid  is  removed.  Although  the  number  of  nodes  in  the  tree  vary,  we  would 
like  each  node  to  contain  a  fixed  amount  of  information  when  representing  a  grid.  Since  the 
number  of  offspring  of  each  node  is  variable,  the  tree  is  implemented  with  each  node  having 
an  offspring  pointer  only  for  the  first  subgrid,  with  sibling  pointers  for  the  rest  of  the  offspr- 
ing. 

(1)  grid  corners 

(2)  number  of  grid  points 

(3)  Ax,  Ay,  At 

(4)  level  in  tree 

(5)  pointer  of  first  offspring  (grid  number) 

(6)  pointer  to  next  grid  at  the  same  level  (grid  number) 

(7)  time  to  which  this  grid  has  been  integrated 

(8)  index  into  main  storage  array  where  approximate  solution  values 
are  stored  at  the  most  recent  time. 

(9)  index  into  main  storage  array  where  approximate  solution  values 
are  stored  at  the  previous  time 

(10)  index  into  main  storage  array  where  perimeter  fluxes 
are  stored  for  conservative  interfaces. 

The  other  data  structure  used  in  our  mesh  refinement  program  manages  the  large  array 
where  the  solution  values  on  all  grids  are  stored.  This  storage  area  is  managed  by  keeping  a 
list  of  used  and  available  blocks  of  storage.  When  a  grid  is  created,  contiguous  storage  space 
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is  reserved  from  the  sorted  list  of  free  blocks  using  a  first-fit  algorithm  [Knuth|.  In  this  algo- 
rithm, the  list  of  free  blocks  is  scanned  until  a  large  enough  block  is  found.  The  requested 
space  is  allocated,  and  the  unused  portion  is  returned  to  the  list.  Reclaimed  space,  which 
occurs  when  a  grid  is  no  longer  needed,  is  inserted  back  into  the  linked  list  of  free  space. 
For  quick  memory  access,  space  is  never  allocated  in  a  circular  fashion  across  the  array  boun- 
dary; all  memory  allocation  must  be  contiguous. 

6.    Interface  Conditions 

Since  each  grid  is  separately  defined,  with  its  own  solution  vector,  a  grid  can  be 
advanced  independently  of  other  grids,  except  for  the  determination  of  its  boundary  values. 
The  boundaries  of  the  fine  grids  will  in  general  be  interior  to  the  problem  domain.  These 
boundary  values  will  usually  come  from  the  coarse  grid.  We  will  not  discuss  here  the  prob- 
lem of  the  boundary  conditions  needed  at  the  boundary  of  the  computational  domain.  Our 
aim  is  to  produce  a  method  to  automatically  handle  the  extra  complexity  that  mesh  refine- 
ment introduces. 

There  are  many  ways  to  get  boundary  values  for  a  fine  grid.  The  different  techniques 
depend  on  whether  the  dependent  variables  are  grid  point  centered  or  cell-centered.  In  addi- 
tion, the  difference  scheme  on  a  coarse  grid  must  sometimes  be  modified  at  the  cell  adjacent 
to  a  fine  grid  if  the  interface  is  to  be  conservative.  New  techniques  are  still  under  develop- 
ment and  being  tested.  A  more  complete  discussion  of  conservation  is  in  [Berger(1987)]. 
Here,  we  will  describe  some  of  the  algorithms  currently  used  at  the  interface  between  a  fine 
and  coarse  grid.  We  use  either  the  Coarse  Mesh  Approximation  Method  (CMAM) 
[Ciment(1971)],  or  interpolation  in  space  and  time  to  get  the  boundary  values. 

We  illustrate  the  CMAM  algorithm  for  fine  grid  boundary  conditions  in  one  space 
dimension  and  time.  In  Figure  6.1,  let  the  coarse  grid  be  on  the  left  of  the  interface,  with 
approximate  solution  values  labelled  v0,v  .,  ,v  -;,...,  as  you  move  left.  Let  the  fine  grid  be 
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V       =  U 
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Fig.  6.1  shows  an  interface  in  one  space  dimension  and  time,  with  grid  point  centered 
variables.  A  stencil  for  the  Coarse  Mesh  Approximation  Method  is  illustrated. 


on  the  right,  refined  by  a  factor  of  2,  with  approximate  solution  values  u0,U|,u:,  .  .  .  ,  as 
you  move  right.  Let  Q  be  an  explicit  finite  difference  approximation  applied  on  the  coarse 
grid, 

The  CMAM  formula  determines  the  missing  boundary  values,  ufl,  0  <  is  2,  using  the  same 
difference  approximation  Q  applied  at  v0  but  with  a  smaller  time  step,  and  with  the  values 
needed  to  the  right  of  v0  supplied  by  u;,«.,  etc.  Note  that  if  the  fine  mesh  finite  difference 
scheme  has  a  stencil  with  more  than  three  points,  additional  boundary  values  will  be  needed. 
If  Q  is  a  second  order  accurate  two  time  level  formula,  the  boundary  values  will  be  formally 
as  accurate  as  the  interior,  coarse  grid  calculation.  A  disadvantage  of  the  CMAM  is  that  it 
requires  the  fine  and  coarse  grids  to  meet  exactly  at  a  grid  point,  unless  spatial  interpolation 
is  used.  However,  in  this  case,  we  might  as  well  interpolate  for  the  boundary  value  itself. 
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Interpolation  is  another  method  for  finding  fine  grid  boundary  values.  The  simplest 
procedure  uses  linear  interpolation  in  time  and  bilinear  interpolation  in  space  for  a  two 
dimensional  computation.  Since  the  coarse  grid  is  integrated  in  time  before  the  fine  grid, 
boundary  values  can  be  obtained  for  all  intermediate  fine  grid  time  steps  by  interpolating  in 
space  and  time  using  only  the  enclosing  coarse  grid  points.  For  example,  in  Fig.  6.1  u",']  : 
would  be  obtained  by  linear  interpolation  using  v'^  and  vq'1.  We  use  this  technique  with 
cell-centered  variables  or  rotated  grids,  since  the  grid  points  of  the  fine  and  coarse  grid  do 
not  coincide. 

Finally,  it  is  sometimes  possible  to  obtain  a  boundary  value  for  a  fine  grid  from  the 
interior  of  an  adjacent  fine  grid. 

6.1.    Conservation 

For  problems  with  shocks,  it  is  desirable  to  use  a  conservative  difference  scheme  at  the 
interface.  We  assume  the  pde  is  in  conservation  form, 

",  -/(«),  =  0 

U(x,t  =  0)    =    Koto 

in  one  dimension,  and  a  conservative  explicit  finite  difference  scheme  of  the  form 

U'J       ~    U'.J  7A f'-l2,J      ti-\2.j) 

is  used  in  the  interior  of  each  grid.  These  considerations  all  generalize  to  two  or  three  spatial 
dimensions.  [Lax  and  Wendroff]  showed  that  schemes  in  conservation  form  converge  to  a 
weak  solution, 

Jfu4>,  -  /(««)<|>t  dx  dt  +  Ju0(x)<t»(x)  dx  =  0, 

for  any  smooth  test  function  <t>(x,r)  with  compact  support.  This  guarantees  that  captured 
shocks  satisfy  the  Rankine  Hugoniot  conditions 
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*M     [fl 

where  the  brackets  denote  the  jump  in  the  quantity  at  the  discontinuity,  and  s  is  the  shock 
speed.  We  sketch  the  Lax-Wendroff  argument  here,  since  we  must  generalize  it  for  the  case 
of  conservation  on  irregular  grids.  Take  the  conservation  form  difference  scheme,  multiply 
by  h,k  and  a  test  function  <J>"  =  4>(jh,nk).  and  sum  over  ail  grid  points  j  and  n  s  0,  giving 

u"~'  -u" 

22  *-r-L*;"> 


j  n&0  * 


j  i»0 


2  j     :  '  "   "  •  "   " — v  '  f"  >  u'  <t>;  kh. 


Applying  summation  by  parts  to  the  left  side  of  this  equation  gives 


1 2  -•'  fc  "•'♦;**  =  2  2  ";v^\    ""m  -  I"?*?''. 

;  <i»0 


which  converges  to  Jj-u<b,dx  dt  -  Ju0(x)<&(x)dx  under  suitable  hypotheses  as  h,k~0.  A 
similar  summation  is  done  on  the  right  hand  side. 

Another  way  to  look  at  the  above  is  that  the  integral  in  the  definition  of  a  weak  solution 
is  approximated  by  the  trapezoid  rule  on  a  uniform  grid  with  mesh  width  h, 

Ju*  dx  =  2";  <f>,  *  +  0(hr). 
j 

The  pde  conserves  the  quantity  /(r)  =  fu(x,i)dx,  except  for  the  flux  at  the  boundaries  of  the 
computational  domain.  The  discrete  approximation  to  /(f), 

S"  =    2  h  "; 

is  conserved  exactly  by  a  numerical  scheme  in  conservation  form.  It  is  easy  to  show  using 
summation  by  parts  that  if  a  conservative  numerical  scheme  is  used  at  all  but  one  point,  and 
the  difference  scheme  at  that  point  is  derived  so  that  the  discrete  integral  S  is  exactly  con- 
served, then  a  weak  solution  is  obtained  [Berger(1987)].  It  is  this  integral  approximation  and 
this  numerically  conserved  quantity  S    that  we  generalize  when  considering  conservation  on 
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irregular  grids,  where  the  mesh  is  refined  in  space  and  time. 

We  illustrate  this  by  making  a  conservative  generalization  of  the  Coarse  Mesh  Approxi- 
mation Method  first.  Let  the  trapezoid  rule  be  used  to  approximate  an  integral  over  a  mesh 
with  an  interface  at  x  =  0.  For  simplicity  of  exposition,  let  the  mesh  be  refined  by  a  factor  of 
2,  so  that 

J  f(x,t)dx  ~     V    f(x,t)2h  +    %f(x,t)h  -  ^/(0,0- 

A  conservative  difference  scheme  is  used  at  all  grid  points  away  from  the  interface.  A  con- 
servative difference  approximation  for  v(0,r)  and  u(0,t),  the  fine  and  coarse  grid  interface 
points,  is  derived  by  insisting  that  the  solution  exactly  conserve  this  discrete  integral  approxi- 
mation at  each  time  step.  The  summation  by  parts  will  give  0(1)  terms  at  the  interface  that 
do  not  in  general  cancel,  since  the  mesh  widths  abruptly  change  there.  However,  the  condi- 
tion that  they  cancel  at  least  to  0(h)  determines  the  conservative  interface  equations,  such  as 

Vq         -  vq         _    f|  :        F  -i  2 
k  3h_ 

2 

,."-1       ..n-12  r"-l  I    _    r" 

vq      ~  vq  n:  ^  -  I  2 

k  }h_ 

2 

1*0  -    v0       . 

which  have  the  interesting  non-symmetric  stencil  shown  for  the  case  of  mesh  refinement  by  a 
factor  of  3  in  Figure  6.2.  One  pleasant  feature  of  this  scheme  is  that  you  don't  have  to  save 
all  the  intermediate  fine  grid  values  in  time,  such  as  uj  " '  2. 

Unfortunately,  for  cell-centered  variables,  we  have  not  yet  found  such  a  simple  conser- 
vative interface  equation.  This  is  because  the  fine  and  coarse  grid  points  are  not  identical, 
even  when  the  grids  are  not  rotated.  Our  procedure  is  to  use  interpolation  boundary  condi- 
tions for  the  fine  grid.    This  determines  a  flux  across  the  interface  from  the  coarse  to  the  fine 
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Fig.  6.2  Conservative  interface  stencil  lor  mesh  refinement  in  time  and  space. 


grid.  We  use  this  to  redefine  the  coarse  grid  flux  at  the  coarse  grid  points  adjacent  to  the 
fine  grid  boundary.  For  example,  in  rwo  dimensions,  the  difference  scheme  at  such  a  coarse 
grid  point  ij  immediately  to  the  right  of  the  fine  grid  should  be 


iU-if 


r-l r-l 


u,j(t+  &tc„ne)  =  «,.,(/)  +  — [f,-i;j(0 rS  2ft-i: 


m  -  p 


(t^q\t,me)\ 


Ay 


[G,,-i:(f)-G;,-i:(r)], 


where  Ax  and  Ay  refer  to  the  coarse  grid  mesh  widths.  The  double  sum  is  due  to  the  refine- 
ment in  time:  for  a  refinement  ratio  of  r,  there  are  r  times  as  many  steps  taken  on  the  fine 
grid  as  the  coarse  grid.  If  the  cell  to  the  north  of  (i,j)  were  also  refined  the  flux  G,  ..  |  ; 
would  be  replaced  by  the  sum  of  fine  fluxes  as  well.  Unfortunately,  to  implement  this  con- 
servative modification,  the  perimeter  fluxes  for  all  fine  grids,  and  the  coarse  grid  fluxes 
corresponding  to  outer  boundaries  of  fine  grids,  must  be  temporarily  saved  each  time  step. 
The  modification  is  a  negligible  amount  of  work,  however,  since  it  only  takes  about  .3%  of 
the  run  time  of  a  typical  run.  This  conservative  fixup  is  what  was  referred  to  in  describing  the 
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updating  step  of  the  adaptive  mesh  refinement  algorithm  in  section  2. 

6.2.    Stability 

In  the  interface  conditions,  as  with  other  non-physical  boundary  conditions,  subtle  non- 
physical  instabilities  can  arise.  These  instabilities  do  not  necessarily  cause  the  solution  to 
blow  up,  but  can  lead  to  persistent  oscillations  near  interfaces  and  slow  convergence  to  steady 
state.  We  will  give  some  examples  of  unstable  interface  conditions  (whose  instability  we 
discovered  the  hard  way,  by  trying  them),  and  list  some  other  interface  conditions  that  we 
can  prove  are  stable.  The  analysis  uses  the  general  stability  theory  of  [Gustafsson,  Kreiss, 
and  Sundstrom].  We  note  that  we  have  always  been  able  to  find  stable  interface  conditions 
with  local  second  order  accuracy  (after  some  work). 


fl^ 


o 


<> 
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f,    =    f0    =    f  /    2 

1  2  coarse   ' 


Fig.  6.3  shows  an  unstable  interface  procedure. 


One  example  of  an  unstable  interface  condition  arose  from  our  steady  flow  computa- 
tions as  follows.  Instead  of  the  conservative  interface  procedure  described  above  for  cell- 
centered  variables,  where  the  coarse  grid  difference  equation  must  be  explicitly  modified, 
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there  is  a  simpler  procedure  one  might  think  of  trying  first.  Compute  the  coarse  grid  flux. 
and  then  divide  it  among  the  fine  grid  cells  bordering  it,  as  illustrated  in  Figure  6.3.  I  his 
leads  to  extremely  slow  convergence,  which  we  eventually  attributed  to  a  mild  radiative  insta- 
bility. Dividing  the  coarse  flux  sets  off  an  odd/even  oscillation  on  the  fine  grid,  which  is  sup- 
ported by  non-dissipative  difference  schemes  like  central  differencing.  Another  instability 
comes  from  using  the  Leapfrog  scheme  with  linear  interpolation  for  the  fine  grid  boundary 
values.  This  was  found  by  Oliger.  We  have  been  able  to  show  that  this  is  unstable  only  for  a 
mesh  refinement  factor  of  2. 

On  the  positive  side,  we  have  shown  analytically  that  Lax-Wendroff  with  linear  interpo- 
lation for  the  fine  grid  boundary  values  is  always  stable,  for  any  mesh  refinement  factor  r 
[Berger(1985)).  We  can  also  prove  that  the  conservative  procedure  described  above  is  stable 
for  mesh  refinement  in  space  only,  with  Runge  Kutta  time  stepping.  Of  course  in  practise  it 
looks  stable  for  mesh  refinement  in  space  and  time,  but  this  is  a  much  more  difficult  theorem 
to  prove. 

7.    Applications 

The  time-dependent  adaptive  mesh  refinement  algorithm  (AMR)  just  described  has 
been  combined  with  FL052  [Jameson)  to  compute  steady  2D  transonic  flow  around  an  air- 
foil. Several  simplifications  of  AMR  are  possible  here  since  we  are  not  interested  in  a  time 
accurate  solution.  A  new  complexity  is  introduced  by  using  multigrid  to  accelerate  the  con- 
vergence of  the  solution  to  steady  state,  since  the  fine  grids  are  adaptively  created  during  the 
solution  process,  and  are  not  global.  We  will  briefly  review  the  well-known  FL052  algo- 
rithm. We  then  describe  the  modifications  necessary  to  interface  AMR  and  FL052.  and  show 
results  of  solving  the  Euler  equations. 

The  full  Euler  equations  in  2D  are  written  in  integral  form. 


j-JJwdxdy  -  J    fdy-gdx  =  0, 
dt      n  all 
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and  ft  is  any  control  volume.  A  finite  volume  discretization  with  cell-centered  variables  is 
used,  where  the  flux  is  evaluated  at  the  boundary  of  a  given  cell  using  the  average  of  the 
values  in  the  adjacent  cells.  On  a  uniform  Cartesian  grid  this  is  the  same  as  central  differenc- 
ing. This  spatial  discretization  procedure  leads  to  a  system  of  ordinary  differential  equations 
of  the  form 


dt 


(hw)  -  Qw-Dw  =  0, 


where  Qw  is  the  approximation  to  the  Euler  terms,  Dw  is  an  added  dissipative  term,  and  h  is 
the  cell  area.  The  dissipation  is  introduced  by  a  combination  of  second  and  fourth  order 
differences,  which  are  switched  on  by  pressure  gradients.  The  same  dissipation  formulas  are 
used  in  the  integration  step  on  each  grid,  except  at  the  boundaries  of  the  fine  grids,  where 
the  fourth  order  stencil  is  too  large.    In  this  case  only  the  second-order  dissipation  is  used. 

The  ODE's  are  integrated  using  a  modified  four  stage  Runge  Kutta  scheme  in  which  the 
dissipative  terms  are  only  evaluated  once.  At  each  time  step  the  solution  is  updated  by  the 
following  sequence: 

ww  =  w(0)  _   *LtQWM-DwM) 
4/:  ' 


w(2)    =   WW   _    £L(QwU)-DwM) 
3/i 


,(3)    =    „(0)    _    Ai(evvC)_DH,'0)) 

2/.  ' 


Ww'     =     H>' 


,(4)    =    w(0)    _    Al(evv(3)_Dvv(0))) 


»VV         =    Wl 


where  w'0)  is  the  value  at  the  beginning  of  the  time  step,  and  w'4)  is  the  final  updated  value. 
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The  time  step  limit  for  this  scheme  is  almost  the  same  as  that  of  a  standard  fourth  order 
Runge  Kutta  scheme.  Other  multistage  Runge-Kutta  schemes  have  been  used  instead,  with 
coefficients  chosen  to  maximize  other  properties  such  as  dissipation  [Turkel),  rather  than 
time  step. 

We  plug  this  well-known  scheme  in  to  the  AMR  algorithm.  The  exterior  of  the  airfoil 
is  discretized  using  an  O  mesh.  In  the  computational  plane,  this  is  simply  a  rectangle.  In 
addition,  the  fine  grids  will  be  in  the  same  coordinate  system  as  the  coarse  grid,  so  that  they 
are  also  rectangles  in  the  computational  plane.  For  example,  the  refinement  at  the  leading 
edge  in  Figure  7.]a  is  the  center  rectangle  in  the  computational  domain  shown  in  Figure  7.1b. 
Since  the  grid  is  periodic  in  the  4  direction,  with  the  break  at  the  trailing  edge,  the  trailing 
edge  refined  grids  are  the  left-  and  rightmost  rectangles  in  Figure  7.1b. 
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Fig.  7.1  Fine  grids  at  the  leading  and  trailing  edges. 


Pseudo-time-stepping  is  used  to  integrate  the  solution,  with  a  variable  At  and  a  fixed 
Courant  number  at  every  grid  point  on  every  grid.  Since  this  is  no  longer  time  accurate,  there 
is  no  reason  to  take  several  smaller  time  steps  on  the  fine  grids  for  every  one  coarse  grid 
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step.  We  simply  take  one  iteration  on  each  grid.  In  addition,  there  is  no  need  to  estimate  the 
error  every  few  time  steps.  Instead,  the  solution  procedure  starts  by  time  stepping  on  a  single 
global  grid.  Since  the  initial  conditions  are  typically  uniform  flow,  we  wait  until  the  solution 
has  settled  down,  typically  with  a  residual  of  ^lO"2,  before  applying  the  error  estimator  and 
adaptively  adding  fine  grids.  The  new  steady  state  solution  is  calculated,  and  again  the  error 
estimator  and  grid  generator  are  invoked  to  see  if  even  finer  grids  are  needed,  or  the  existing 
fine  grids  need  to  be  re-located.  This  adjustment  is  necessary  because  the  resolution  at  the 
leading  edge,  for  example,  can  sufficiently  change  the  shock  location  so  that  it  moves  out  of  a 
refined  grid  region.  We  continue  using  the  Richardson  error  estimates  to  determine  how  fine 
a  grid  is  necessary  to  obtain  a  specified  accuracy  in  the  solution. 

Finally,  the  error  estimation  procedure  itself  can  be  simplified,  since  only  the  spatial 
discretization  error  needs  to  be  measured.  Once  the  solution  is  nearly  converged,  one  resi- 
dual calculation  using  every  other  grid  point  in  the  mesh  gives  is  the  local  truncation  error  in 
space.    This  is  similar  to  the  tau  extrapolation  procedures  used  in  multigrid  computations. 

Things  become  more  complicated  with  the  introduction  of  the  multigrid  acceleration 
strategy.  In  this  scheme,  the  equations  on  a  coarse  grid  are  augmented  by  a  forcing  function. 
If  wh  represents  the  solution  on  a  fine  grid,  and  w2h  the  solution  on  a  coarse  grid,  then  the 
coarse  grid  equations  are 

w$  =  w#>  -  Ataw  (  Rw'{l<  -  Plh) 

wfl  =  w<$  -  Afa<h  (Rw'$  -  />:„) 

etc,  where  Pzh  =  ^R(*>h)~R(w*fti)-  The  sum  is  taken  over  the  four  fine  cells  comprising 
one  coarse  cell.  This  formulation  insures  that  when  the  residual  is  zero  on  the  fine  grid,  the 
solution  on  the  coarse  grid  does  not  change.  After  one  iteration  on  a  coarse  grid,  the  correc- 
tion is  interpolated  back  to  the  fine  grid  using  bilinear  interpolation. 
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Because  of  the  generality  of  the  data  structures  in  AMR.  very  little  has  to  change  to 
include  multigrid.  The  updating  procedure  of  AMR  is  the  same  as  the  restriction  operator  in 
multigrid.  However,  now  both  the  solution  and  the  residual  are  transferred  from  the  fine  to 
the  coarse  grid.  To  use  multigrid  on  fine  grids  that  are  not  global,  means  that  some  coarse 
cells  have  a  forcing  function  from  a  finer  grid  and  some  do  not.  This  is  easily  arranged, 
since  the  AMR  strategy  previously  described  already  accounts  for  the  fact  that  some  coarse 
solution  values  are  updated  form  a  finer  grid  and  some  are  not.  The  only  new  operation  is 
prolongation:  the  interpolation  of  the  correction  to  the  coarse  grid  variables  back  onto  the 
fine  grid. 

The  use  of  multigrid  on  patched  grids  however  is  tricky.  The  implementation  of  the 
boundary  conditions  for  the  fine  grids  and  conservation  modification  of  the  coarse  grid 
affects  the  convergence  rate  of  the  solution  to  steady  state.  We  currently  use  the  following 
procedure,  but  caution  that  this  work  is  still  in  progress. 

A  four  stage  Runge  Kutta  scheme  used  with  central  differencing  in  space  needs  one 
boundary  value  at  each  stage,  or  four  boundary  values  at  the  initial  stage.  We  use  the  former 
alternative,  where  the  boundary  value  is  obtained  from  the  coarse  grid.  This  value  is  held 
fixed  for  all  stages,  since  at  steady  state,  there  is  no  change  in  the  solution  at  each  stage. 
However,  this  simplification  leads  to  slower  convergence  than  is  observed  when  the  fine  grids 
are  global. 

On  the  coarse  grid,  some  difficulty  arises  at  grid  points  adjacent  to  fine  grid  boundaries. 
We  give  a  simple  example  illustrated  in  Figure  7.2.  At  steady  state,  the  conservative  differ- 
ence scheme  at  the  point  u,  ,  on  the  coarse  grid  uses  a  stencil  based  on  values  at 
"'-i.;'"'j-i'M'j-i.  and  the  points  v2 .*  and  v:<-i  from  the  fine  grid.  During  a  coarse  grid 
iteration,  v:t  and  v;t.,  don't  change.  In  addition,  any  change  in  the  coarse  grid  variables 
"underneath"  the  fine  grid  won't  be  felt  until  the  next  coarse  grid  iteration.  We  call  this  a 
"no-flow"  interface  equation. 
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Fig.  7.2  illustrates  the  interfaces  conditions  between  the  fine  and  coarse  grid. 


This  situation  can  be  remedied  by  using  the  following  "flow-through''  equations  at  these 
coarse  grid  points.  We  implement  the  conservation  modification  via  the  forcing  function,  so 
that  the  difference  scheme  at  u,j  does  use  the  point  u,-\j.  Define  the  forcing  function  spe- 
cially for  coarse  points  adjacent  to  the  fine  grid,  so  that  (illustrating  with  a  linear  wave  equa- 
tion) Pu  =  -", -u  +  (v:.t  +  v2.v-i)/2-  At  steady  state,  conservation  is  still  maintained,  but 
on  multi-stage  Runge  Kutta  iterations,  after  the  first  stage,  the  coarse  grid  points 
u,-\j,u,-2j,  etc.  can  play  a  role  if  they  have  changed.  This  version  of  the  AMR  multignd 
uses  approximation  10%  fewer  iterations  than  the  non-flow  through  boundary  conditions 
above. 

Further  evidence  that  the  interface  conditions  are  still  slowing  convergence  is  in  the  fact 
that  three  stage  Runge  Kutta  schemes  suffer  less  of  a  slowdown  than  five  stage  schemes,  for 
example.  We  are  continuing  to  work  on  better  equations  for  the  interface  conditions.  How- 
ever, even  with  the  loss  in  convergence  rate,  the  AMR  multigrid  still  yields  computational 
savings  over  global  grid  multigrid  schemes,  owing  to  the  smaller  fine  grid  sizes  and  reduced 
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cost  per  fine  grid  iteration. 

The  following  figures  compare  the  results  of  computing  the  solution  on  a  uniform 
coarse  grid,  a  coarse  grid  with  fine  grid  patches,  and  a  globally  refined  grid.  In  all  cases,  we 
obtain  the  same  accuracy  as  the  solution  on  a  uniform  fine  grid  by  refining  only  a  fraction  of 
the  coarse  grid.  The  figures  are  taken  from  [Berger  and  Jameson]. 

The  first  test  case  is  for  non-lifting  subsonic  flow  over  a  NACA  0012  airfoil.  The  Mach 
number  is  .500  with  zero  degrees  angle  of  attack.  Figure  7.3  shows  the  pressure  coefficient 
calculated  on  a  32  by  8  coarse  grid.  The  drag  coefficient  is  .0049.  Figure  7.4  shows  the 
pressure  coefficient  calculated  on  a  64  by  16  grid.  The  drag  coefficient  is  .0011.  Figure  7.5 
shows  the  solution  using  a  grid  of  size  128  by  32  grid,  with  a  drag  coefficient  of  .0002.  The 
drag  coefficient  is  converging  like  h2  to  its  expected  value  of  zero.  Figure  7.6  shows  a 
refined  grid  solution  based  on  a  32  by  8  underlying  coarse  grid,  with  refined  grid  patches 
shown  in  Figure  7.7.  The  drag  coefficient  in  this  case  is  .0009.  Figure  7.8  shows  a  refined 
grid  solution  based  on  a  64  by  16  underlying  coarse  grid.  The  refined  grid  for  this  case  is 
shown  in  figure  7.9.  The  drag  coefficient  is  reduced  to  .0001.  In  both  cases,  the  accuracy  of 
the  solution  on  the  uniformly  next  finer  level  grid  is  recovered  by  using  small  grid  patches  at 
the  leading  and  trailing  edges  of  the  airfoil.  More  dramatic  savings  can  be  achieved  by  using 
three  nested  levels  of  grid  refinement,  as  shown  in  figure  7.10.  The  solution,  with  a  drag 
coefficient  of  .0003,  is  shown  in  figure  7.11. 

The  second  test  case  is  transonic  flow  containing  a  shock  wave.  Figure  7.12  shows  the 
pressure  coefficient  for  a  NACA  0012  airfoil  at  Mach  .8  with  zero  degrees  angle  of  attack. 
The  mesh  used  for  this  computation  is  64  by  16  cells.  When  the  grid  is  refined  (using  an 
error  tolerance  of  .005),  as  shown  in  figure  7.13,  the  solution  obtained  is  almost  identical  to 
the  solution  computed  on  a  128  by  32  mesh  (compare  figures  7.14  and  7.15).  In  the  coarse 
grid  run,  the  entropy  behind  the  shock  was  computed  to  be  .0072,  in  the  multiple  grid  run  it 
was  .0052,  and  in  the  fine  grid  run  the  entropy  was  .0054.    In  this  mesh  refined  solution, 
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21%  of  the  coarse  grid  was  refined  by  a  factor  of  2  in  both  coordinate  directions.  If  the  error 
tolerance  for  mesh  refinement  is  less  stringent  (.025),  so  that  only  the  leading  and  trailing 
edges  are  refined  (as  in  figure  7.16),  the  solution  is  only  slightly  worse  (figure  7.17).  The 
entropy  production  across  the  shock  is  .0046  in  this  case.  In  this  run  only  10%  of  the  coarse 
grid  is  refined.  This  also  shows  that  most  of  the  error  in  computing  the  drag  is  due  to  the 
lack  of  resolution  at  the  leading  and  trailing  edges. 

The  final  computations  we  show  come  from  a  study  of  shock  reflection  off  a  wedge. 
The  incident  shock  is  travelling  at  mach  10  up  a  ramp  at  an  angle  of  30  degrees,  and  y  =  1.4 
for  the  gas.  These  computations  use  the  time  dependent  version  of  AMR,  with  mesh  refine- 
ment by  a  factor  of  4.  The  figures  are  taken  from  [Berger  and  Colella].  Figure  7.18  shows  a 
contour  plot  of  density.  The  contour  plotter  uses  solution  values  taken  from  the  finest  avail- 
able grid.  One  can  tell  by  looking  at  the  reflected  shock  where  we  have  turned  off  the  grid 
refinement  -  the  shock  width  of  the  reflected  shock  increases  by  an  order  of  magnitude.  Fig- 
ure 7.19  shows  the  fine  grids  around  a  blowup  of  the  mach  stem.  Three  levels  of  grid  refine- 
ment were  used.  80%  of  the  run  time  was  spent  in  the  integrator;  95%  of  this  integration 
time  was  spent  on  the  finest  level  grids.  Since  these  grids  cover  only  10%  of  the  domain,  the 
computational  savings  in  using  adaptive  mesh  refinement  is  roughly  a  factor  of  8.  The  coarse 
grid  in  this  run  is  100  by  20.  A  uniformly  refined  grid  run  at  the  same  level  of  refinement  as 
the  finest  grids  would  have  used  2  million  words  of  memory.  The  adaptively  refined  compu- 
tation used  a  maximum  of  830,000  words  during  the  run.  The  savings  in  storage  are  less  than 
the  savings  in  CPU  time  since  our  present  code  is  rather  profligate  with  storage.  We  feel 
that  on  current  computers  it  is  more  important  to  reduce  the  CPU  time  than  the  memory 
requirements. 
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8.    Conclusions 

Let  us  summarize  by  discussing  some  of  the  weaknesses  in  the  AMR  approach.  It  can 
be  difficulty  to  find  conservative  accurate  boundary  conditions  that  minimize  wave  reflection 
at  grid  interfaces.  However,  this  has  only  been  a  problem  with  non-dissipative  difference 
schemes.    We  also  have  no  experience,  as  yet,  with  using  implicit  difference  schemes. 

At  present,  we  align  fine  grids  with  the  coarse  grid  only  because  we  have  not  yet 
worked  out  conservative  boundary  conditions  for  use  with  rotated  fine  grid  patches.  We 
would  be  happier  if  we  could  find  simple  and  accurate  conservative  interface  conditions  for 
rotated  grids.  Rotation  would  have  the  advantage  of  better  aligning  the  fine  grid  with 
features  of  the  flow  (e.g.  shocks),  and  minimizing  the  area  of  the  coarse  grid  that  needs  to  be 
refined.  In  the  future,  however,  if  we  combine  local  grid  refinement  with  moving  mesh 
methods,  the  use  of  rotated  rectangles  may  become  less  desirable.  [Gropp]  has  started  in  this 
direction  with  his  work  using  embedded  grids  that  are  aligned  with  the  coarse  grid,  but  allow- 
ing the  fine  grids  to  move  in  any  direction.  We  are  also  beginning  to  combine  AMR  with  the 
Brackbill  and  Saltzman  variational  approach  to  mesh  movement.  We  mention  that  AMR  has 
already  been  combined  with  the  shock  fitting  algorithm  of  [Chern  and  Colella]. 

The  complexity  of  the  AMR  code  might  be  intimidating  to  a  new  user.  Our  3000  lines 
of  Fortran  must  be  contrasted  with  the  1600  lines  of  FL052.  for  example.  However,  a  big 
code  is  not  necessarily  a  fragile  code.  We  have  been  careful  to  develop  AMR  to  make  it 
automatic  and  robust.  In  addition,  a  user  should  be  able  to  use  AMR  without  having  to 
understand  it  all.  This  is  the  way  we  use  compilers,  for  example.  This  makes  it  important  to 
develop  AMR  in  a  modular  way.  A  user  should  be  able  to  plug  in  an  integrator  for  a  new 
problem  without  knowing  details  about  how  the  more  computer  science  oriented  parts  of 
AMR  work,  but  knowing  that  these  other  parts  will  indeed  work.  We  have  already  demon- 
strated this  modularity  by  using  AMR  in  conjunction  with  FLOS2,  with  PPM  [Colella  and 
Woodward;  Colella],  and  we  are  currently  working  on  a  problem  in  magneto-hydrodynamics. 


-    40   - 


CP 


MACH  0.500ALPHA  0. 

CL   O.OOOOCD   0.0049CM  -0.0000 


* 
* 


Fig."Vi  Pressure  coefficient  using  a  32  x  8  coarse  grid. 


CP 


MACH  0.500ALPHA  0. 

CL      O.OOOOCD      O.OOUCM    -C .  OOOC 

********* 


* 


Kip  1  S    Pressure  coefficient  using  a  W  x  16  coarse  grid. 


CP 


MACH  0.500ALPHA  0. 

CL   O.OOOOCD   0.0002CM  -0.0000 


* 
* 


Fig. "IS   Pressure  coefficient  using  a  128x32  coarse  grid. 
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Fig. It    Pressure  coefficient  using  the  grid  in  Fig.  9. 
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Fig.  >}    Pressure  coefficient  using  the  grid  in  Fig.  11. 
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Fig.^-.lC   32  x  8  coarse  grid  with  two  levels  of  refinement. 
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Fig.1l?    Pressure  coefficient  using  a  64  x  16  coarse  grid. 
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Fig  Vila    64  x  16  coarse  grid  with  one  level  of  refinement,  using  a  less 
slringenl  error  tolerance. 
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Fig.  7.18  shows  a  contour  plot  of  the  density. 


Fig.  7.19  shows  a  blowup  of  the  region  containing  fine  grids  at  level  2  and  3. 
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